Anomalous melting behavior of solid hydrogen at high pressures 
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Hydrogen is the most abundant element in the universe, and its properties under conditions of 
high temperature and pressure are crucial to understand the interior of large gaseous planets and 
other astrophysical bodies. At ultra-high pressures solid hydrogen has been predicted to transform 
into a quantum fluid, because of its high zero-point motion. Here we report first-principles two- 
phase coexistence and Z-method determinations of the melting line of solid hydrogen in a pressure 
range spanning from 30 to 600 GPa. Our results suggest that the melting line of solid hydrogen, as 
derived from classical molecular dynamics simulations, reaches a minimum of 367 K at ~430 GPa; 
at higher pressures the melting line of the atomic Cs-IV phase regains a positive slope. In view of the 
possible importance of quantum effects in hydrogen at such low temperatures, we also determined 
the melting temperature of the atomic Cs-lV phase at pressures of 400, 500 and 600 GPa, employing 
Feynman path integral simulations. These result in a downward shift of the classical melting line by 
~100 K, and hint at a possible secondary maximum in the melting line in the region between 500 
and 600 GPa, testifying to the importance of quantum effects in this system. Combined, our results 
imply that the stability field of the zero-temperature quantum liquid phase, if it exists at all, would 
only occur at higher pressures than previously thought. 

PACS numbers: 62.50.-p, 64.70.dj, 65.20.-w, 65.40.-b 
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The characterization of the physico-chemical proper- 
ties of hydrogen at high pressures and temperatures is of 
fundamental interest in physics, astrophysics, chemistry 
and planetary science. Previous theoretical studies [l|, Q 
have predicted that at sufficiently high pressures solid 
hydrogen would transform into a quantum liquid phase 
that would constitute a new state of matter, due to the 
fact that hydrogen has a very large zero-point energy at 
ultra-high pressures. Results from first-principles molec- 
ular dynamics (MD) simulations have indicated that liq- 
uid hydrogen undergoes a molecular to atomic transfor- 
mation fsj, which is predicted to occur at a pressure of 
125±10 GPa along the 1500 K isotherm. Another study, 
employing similar techniques, has reported the finding 
of reentrant behavior in the melting of hydrogen, i.e. 
the existence of a maximum along the melting line, pre- 
dicted to occur at '^80 GPa and '^900 K, followed by a 
decrease of the melting temperature at higher pressures. 
This prediction has been subsequently confirmed by ex- 
periments 043' which place the maximum at ~106 GPa 
and 1050±60 K. The occurrence of reentrant melting be- 
havior is understood to be a necessary, though by itself 
not sufficient, condition for the existence of the high pres- 
sure quantum liquid phase predicted by Ashcroft [4I , and 
thus it was taken as a strong indication of its existence by 
Bonev and coworkers With the assumption that the 
negative slope of the melting line persists to high enough 
pressures, it was estimated that the quantum fluid 
state would occur at pressures close to 400 GPa. Recent 
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first-principles MD and quantum Monte Carlo [8j simu- 
lations have confirmed once more the reentrant behavior 
of the melting line, and determined a coexistence point 
at which both the molecular and atomic fluids coexist 
with the solid (phase I); this coexistence point could ei- 
ther be a triple point or, if the insulator-metal transition 
occurring in the liquid actually extends into the solid, a 
quadruple point. 

Here we report first principles calculations of the melt- 
ing temperatures of five different crystalline phases of 
solid hydrogen (four molecular ones and one atomic), 
which are the energetically most competitive phases in 
the pressure range 30-600 GPa. These phases are the 
molecular d-hcp (phase I), partially ordered hep {po-hcp 
or phase IV), Cmca-4, Cmca-12 and quasi-molecular 
mC24 phases, plus the atomic Cs-IV phase. We have 
used two reliable methods (two-phase coexistence and 
the Z-method [l3|) to predict the melting curve of hydro- 
gen in this range of pressures. Our results indicate that 
the high-pressure Cs-IV phase of hydrogen has a melt- 
ing line with a positive slope. When considered together 
with the reentrant behavior of the melting line at lower 
pressures, this results in the melting line having a mini- 
mum value of 367 K at a pressure of ~432 GPa. While 
the occurrence of this minimum along the melting line 
does not preclude the existence of the predicted quantum 
liquid phase 0, it would nevertheless push its stability 
field to higher pressures than previously thought. How- 
ever, these results are based on classical simulations, and 
due to the lightness of hydrogen and the relatively low 
melting temperatures that result in this range of pres- 
sures, it is necessary to consider the possible influence 
of quantum effects. This is a point to which we return 
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below. 

The two-phase coexistence approach consists of plac- 
ing within the same simulation box the solid and liq- 
uid phases directly in contact through an interface. The 
melting temperature is determined as that at which both 
phases are seen to remain stably in coexistence. At 
the first-principles level, such simulations are computa- 
tionally demanding, as they require large cells and long 
simulation times, but nevertheless they are accessible 
with modern supercomputer facilities as recently demon- 
strated for the cases of Li JTj and Fe [12] ■ In this work 
we have used first principles molecular dynamics simu- 
lations based on density functional theory (DFT). For 
this task we employed the Vienna ab initio Simulation 
Package (VASP) [13]. Simulations were conducted in the 
NPH ensemble (constant number of particles, N, con- 
stant pressure, P, and constant enthalpy, H) following 
the algorithm of Souza and Martins [l4| as recently im- 
plemented in VASP [III, At a given external pres- 
sure, the temperature adjusts spontaneously and tends 
towards the coexistence temperature, as the unstable 
phase in the simulation is gradually consumed in favor 
of the stable one. By suitably repeating such simulations 
at various temperatures it is possible to obtain the co- 
existence temperature as the temperature at which both 
phases stably coexist and no drift in the instantaneous 
temperature is observed. The presence of both solid and 
liquid phases within the simulation box can be easily 
corroborated either by direct inspection or by plotting 
the particle density calculated at a series of planes par- 
allel to the interface [see Fig ([1])]. Crystal planes ap- 
pear as spikes of high density, whereas no such spikes 
are present in the liquid region. The simulation was 
constrained to retain a tetragonal shape, with the two 
short sides (parallel to the plane of the interface) be- 
gin of equal length. System sizes included 980, 960 and 
960 hydrogen molecules and 2048 hydrogen atoms in the 
(7 X 7 X 10, 2 X 2 X 5, 8 X 8 X 8, 6 X 4 X 10) simulation 
cells of the d-hcp, po-hcp, Cmca-4 and Cs-IV phases, 
respectively. Initial configurations containing approxi- 
mately equal amounts of the solid and liquid phases, 
plus the interface, were then generated. When both 
phases were found to coexist stably for a minimum of 
5 ps they were assumed to be in equilibrium at the tem- 
perature and pressure conditions of the simulation. In 
these simulations the Brillouin zone was sampled with a 
grid of 2 X 2 X 1 fc-points; the self-consistency criterion 
required on the total energy was a variation smaller than 
2 X 10~^ eV between two consecutive iterations, and the 
time step for integration of the equations of motion in 
MD was 0.5 fs. 

The Z method was developed by Belonoshko and col- 
laborators. It was first employed with empirical poten- 
tials (loj . and more recently in combination with first- 
principles methods [1^1 . This method relies on the ob- 
servation that, in the limit of superheating, i.e. when 
the solid is heated to the maximum temperature that 
it can sustain as a metastable phase without transform- 
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FIG. 1: (Color online) The main panel shows the particle 
density of the Cs-IV phase in coexistence with liquid hydro- 
gen at 500 GPa. Regularly spaced density peaks reveal the 
presence of crystal planes on the right, whereas the small den- 
sity oscillations around a constant value seen on the left are 
characteristic of the liquid phase. The upper inset, (a), shows 
an instantaneous configuration resulting from the same coex- 
istence simulation. Crystal planes are clearly visible on the 
right side of the image, and the interface and liquid are dis- 
cernible on the left. The lower inset, (b), shows the time de- 
pendence of the instantaneous temperature of the simulation 
shown in the main panel. The temperature has no tendency 
to drift, oscillating around an average value of 378 K with a 
standard deviation smaller than 20 K, as obtained from the 
last 7 ps of simulation. 



ing into the liquid, its internal energy is equal to that 
of the liquid phase at the thermodynamic melting tem- 
perature. Thus, the problem of determining the melting 
temperature can be rephrased into the the problem of 
determining the temperature of the superheating limit. 
Then, running a microcanonical simulation of the per- 
fect solid at this temperature should result in the liquid 
phase at the equilibrium melting temperature. We have 
applied this technique to calculate the melting line of hy- 
drogen using systems consisting of 48 molecules for the 
molecular phases {po-hcp, Cmca-4, Cmca-12), and 128 
atoms for the atomic Cs-IV phase, employing simulation 
times of at least 10 ps in each case. In these simulations 
the Brillouin zone was sampled a 4 x 4 x 4 fc-point grid, 
and, like for the coexistence simulations, the time step 
used was 0.5 fs. 

Our calculated melting line for the d-hcp structure pro- 
vides a higher melting line than that predicted in previ- 
ous theoretical work j4j in the range 30-200 GPa, but 
nevertheless in good agreement with experimental mea- 
surements @, El- Moreover, we find a melting tempera- 
ture for the po-hcp phase of 562 K at 300 GPa, which is 
similar to that obtained recently [1] (550 K at 290 GPa) 
using free energy calculations combined with quantum 
Monte-Carlo methods. A further test of the reliability 
of our results is the good agreement that we obtain be- 
tween the melting temperatures derived from the coexis- 
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tence calculations and those resulting from the Z-method 
simulations (e.g. the simulated melting temperatures of 
Cs-IV phase at 400 GPa are 350 K and 362 K by using 
Z-method and two-phase coexistence, respectively). 

Recently, phase IV has been found at room tempera- 
ture and pressures above 220 GPa [13, [11] ■ This new 
phase has been interpreted as a mixed structure that 
includes two different molecular layers ^19i] and a par- 
tially ordered hep (po-hcp) structure [1^. Previous the- 
oretical studies [2Q| by the present authors suggest that 
the Cmca-12 structure is less stable than the Cmca-A 
above 250 GPa, as the former has a higher zero-point en- 
ergy than the latter. In view of this, we have calculated 
the melting line of the po-hcp, Cmca-12 and Cmca-A 
phases in the range 200-500 GPa using the ah initio Z- 
method. While determining the melting temperature of 
the Cmca-A and Cmca-12 phases at 300 and 350 GPa, we 
found that these phases transform into the po-hcp struc- 
ture, which indicates that the latter structure is more 
stable than the other molecular structures at tempera- 
tures close to melting in this range of pressures. This re- 
sult also indicates that phase Cmca-12 probably has no 
stability field in the phase diagram of hydrogen, which 
is in good agreement with our previous calculations |20| . 
Note that the melting temperature of the po-hcp (795 K) 
and d-hcp (815 K) structures is very similar at 200 GPa, 
and that the melting line of the po-hcp phase almost co- 
incides with the extrapolation of that of the d-hcp phase. 
This is consistent with the recent experimental observa- 
tion that the d-hcp structure becomes partially ordered, 
transforming into the po-hcp structure at room temper- 
ature and -220 GPa [iTllTsj. 

At higher pressures hydrogen is expected to transform 
into an atomic phase. Our previous studies [2l| have pre- 
dicted the existence of a quasi-molecular mC24 phase in 
the pressure range lying between the stability fields of the 
molecular Cmca-A and the atomic Cs-IV phases. The 
mC24 structure has two different bond lengths (—0.86 
and 0.9 A at 500 GPa) which are longer than that present 
in the molecular Cmca phase (0.78 A at 400 GPa). The 
mC24 phase will dissociate into the atomic phase upon 
temperature increase (at —350 K), due to thermal fluctu- 
ations. Thus, in order to characterize the melting behav- 
ior of hydrogen at pressures in the range 350-600 GPa, 
we have calculated the melting temperatures of the mC24 
and Cs-IV phases by employing both the two-phase co- 
existence and Z-methods. We find that the intersection 
of the melting lines of the molecular and atomic phases 
takes place at -432 GPa and 367 K (see Fig. [^, and 
that the melting temperature of the mC24 structure is 
systematically lower than that of the Cs-IV phase, in- 
dicating that the latter is more stable at temperatures 
close to melting. 

The most intriguing observation to be extracted from 
our melting temperature calculations is that the phase 
diagram of hydrogen should have a second extremum 
along the melting line. The first extremum is the max- 
imum originally reported by Bonev and coworkers 
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FIG. 2: (Color online) Phase diagram of hydrogen. The red 
and violet lines indicate the melting line obtained from our 
two-phase coexistence and Z-method simulations, where T 
and Z mean two-phase coexistence and Z-method, respec- 
tively. The solid symbols are mel ting data as obtained from 
previous experimental work [3. [tI. [23^ [2^ . Theoretical melt- 
ing data from refs. [1, is represented by open symbols. The 
boundaries between phases I-II-III-IV at low temperatures 
are taken from refs. ,25-28]. Inset shows the calculated shear 
modulus of Cmco-4 and Cs-IV structures as a function of 
pressure at 300K. 



and later confirmed by experimental measurements 0-[3l . 
The second extremum we find here occurs at the inter- 
section of the melting lines of the molecular and atomic 
phases (—432 GPa, 367 K, see Fig. [2]) and is a mini- 
mum, rather than a maximum. Although the molecular 
phase has a melting line with a negative slope at pres- 
sures above — 100 GPa, we find the atomic Cs-IV phase 
to have a melting line with a shallow, but nevertheless 
positive slope above 360 GPa. Therefore, the intersection 
of these two melting lines occurs at a minimum. Thus 
it is seen that the melting behaviors of the molecular 
and atomic phases are qualitatively different. In order 
to further assert this difference in melting behavior, we 
have calculated the finite temperature shear modulus of 
the Cmca-A: and Cs-IV phases at 300 K at pressures be- 
tween 300 and 500 GPa. According to the Born melting 
criterion [2^ . melting occurs when the temperature is 
such that the shear modulus of the solid phase reduces to 
zero. Thus we expect to observe a strong correlation be- 
tween the behavior of the shear modulus and that of the 
melting temperature as a function of pressure. The finite 
temperature shear modulus was calculated by performing 
ah initio MD simulations of suitably deformed supercells 
of each phase containing 96 (in the case of the Cmca-A 
structure) or 128 atoms (for the Cs-IV structure). Each 
run consisted of a total of 4000 steps with a time step of 
0.5 fs, and the stress components were averaged over the 
last 2000 steps. In the case of the Cmca-A phase we also 
run the simulation at 400 GPa for a total of 8000 time 
steps, but the averaged stress components were identi- 
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cal to those obtained in the shorter run. The calculated 
shear moduli for these phases at 300 K are shown in the 
inset of Fig. [21 where it can be seen that the moduli of 
the molecular phases decrease with increasing pressure, 
in consonance with the behavior of their corresponding 
melting lines. The shear modulus of the Cs-lV atomic 
phase, in contrast, displays the opposite trend, which is 
again in accordance with its melting behavior. 

Up to this point all our results have been based on clas- 
sical simulations. However, hydrogen being the lightest 
of elements, and the fact that the melting temperatures 
predicted at high pressures are relatively low, make it 
necessary to consider the possible influence of quantum 
effects on the melting temperature. To address this is- 
sue in the particular case of the Cs-IV atomic phase, we 
have employed path integral simulations [29] combined 
with the hysteresis method [30] to determine the melting 
temperature. The hysteresis method, similar in spirit to 
the Z-method discussed earlier, has been shown to be ca- 
pable of producing accurate estimations of the melting 
temperature [13, HH . We employed this method with a 
72-atom super-cell of the Cs-IV structure with 8 beads 
in the quantum ring polymer (doubling the number of 
beads to 16 did not affect the resulting melting temper- 
ature at 400 GPa), and using 5x5x3 fc-points for each 
bead. The resulting melting temperature at 400 GPa is 
^250 K, which is lower than the temperature obtained 
from classical molecular dynamics (350 K), confirming 
the importance of quantum effects in this system. Fur- 
thermore, it can be seen in Fig. ^ that the melting line 
including quantum effects for the Cs-IV is nearly flat in 
the range 400-500 GPa, but seems to regain a negative 
slope above 500 GPa. This could indicate the presence of 
a second local maximum in the melting line of hydrogen, 
followed by the re-establishment of reentrant behavior at 
pressures in the range 500-600 GPa. The existence of a 
negative slope in the melting line in this range of pres- 
sures could again be taken as a possible indication of the 
existence of a quantum liquid phase at high pressure, but 
according to our results we do not expect such a phase, 
if it exists, to be found at pressures below 600 GPa, the 
maximum pressure considered in this study. In this re- 



spect, we note that a very recent study [32| employing 
path integral simulations has reported melting tempera- 
tures of the order of 50 K from 900 GPa onwards. 

In summary, we have calculated the melting line of 
hydrogen up to 600 GPa using first principles molecu- 
lar dynamics methods, considering melting from the d- 
hcp, po-hcp, Cmca-4:, Cmca-12 molecular phases and 
the Cs-IV atomic phase. The obtained classical melt- 
ing line of molecular hydrogen is in good agreement with 
previous experimental and theoretical results. We pre- 
dict, however, that the melting line reaches a minimum 
at a pressure of ~432 GPa, coinciding with the point 
at which the melting lines of the Cmca-A and atomic 
Cs-IV phases cross, and that the classical melting line 
retains a positive (though shallow) slope at least up to 
600 GPa. Quantum effects incorporated approximately 
via Feynman path integral simulations result in a down- 
ward shift of the melting line by as much as '^100 K, and 
hint at the possibility of a secondary maximum along the 
melting line in the region 500 to 600 GPa. These results 
suggest that the proposed low-temperature quantum liq- 
uid phase, if it exists at all, will not be found at pressures 
below 600 GPa. 
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